Recall that the simple linear regression model is:

\[ y_i = \beta_0 + \beta_1 x_i + \epsilon_i \quad \text{for} ~ i=1,\ldots,n, \]

where the error terms \(\epsilon_1,\ldots,\epsilon_n\) are mean zero, independent, and identically distributed random variables.

Let \(\hat{\beta}_1\) denote the parameter estimate of \(\beta_1\). The bias \(\hat{\beta}_1\) in estimating \(\beta_1\) is defined as the expected (i.e., mean) value of \(\hat{\beta}_1 - \beta_1\). A \(100(1-\alpha_1)\)% confidence interval for the slope \(\beta_1\) (assuming independent, normally distributed errors) is \(\hat{\beta}_1 \pm t_{\alpha_1/2,\text{df}} ~ \hat{\sigma}_{\hat{\beta}_1}\), where \(\text{df}\) is the residual (a.k.a., error) degrees of freedom, \(\hat{\sigma}_{\hat{\beta}_1}\) is the standard error of the slope estimate, and \(t_{\alpha_1/2,df}\) is the \(\alpha_1/2\) upper quantile of the \(t\)-distribution with \(\text{df}\) degrees of freedom. The coverage of a confidence interval estimator is the long-run relative frequency that the procedure generates intervals that contain the true value.

The regression coefficients \(\beta_0\) and \(\beta_1\) can be estimated using the method of least-squares. In R, if x and y are numeric vectors of the same length, the least-squares parameter estimates, standard errors, degrees of freedom, etc. are obtained using:

fm <- lm(y ~ x)
fm
summary(fm)

While theory provides results on the bias and coverage in simple linear regression models, the goal is to perform a simulation study to investigate:

The following function conducts such a simulation study. The arguments to the function are:

(x <- seq(-2, 2, length = 50))
 [1] -2.00000000 -1.91836735 -1.83673469 -1.75510204 -1.67346939 -1.59183673 -1.51020408 -1.42857143
 [9] -1.34693878 -1.26530612 -1.18367347 -1.10204082 -1.02040816 -0.93877551 -0.85714286 -0.77551020
[17] -0.69387755 -0.61224490 -0.53061224 -0.44897959 -0.36734694 -0.28571429 -0.20408163 -0.12244898
[25] -0.04081633  0.04081633  0.12244898  0.20408163  0.28571429  0.36734694  0.44897959  0.53061224
[33]  0.61224490  0.69387755  0.77551020  0.85714286  0.93877551  1.02040816  1.10204082  1.18367347
[41]  1.26530612  1.34693878  1.42857143  1.51020408  1.59183673  1.67346939  1.75510204  1.83673469
[49]  1.91836735  2.00000000
intercept <- 57
slope <- 112
true_y <- intercept + slope * x

plot(x, true_y, 'l')

error <- rnorm(length(x), sd = 10)
simulated_y <- true_y + error
plot(x, true_y, type = 'l')
points(x, simulated_y)

(lm_out <- lm(simulated_y ~ x))

Call:
lm(formula = simulated_y ~ x)

Coefficients:
(Intercept)            x  
      59.32       111.89  
(bias <- lm_out$coefficients[2] - slope)
         x 
-0.1086681 
alpha <- .1
ci <- lm_out$coefficients[2] + c(-1, 1) * qt(1 - alpha, df = lm_out$df.residual) * summary(lm_out)$coefficients[2, "Std. Error"]
c(
  lower = ci[1],
  estimate = lm_out$coefficients[2],
  upper = ci[2]
)
     lower estimate.x      upper 
  110.1580   111.8913   113.6246 
library(tidyverse)
library(gganimate)

animation_data <- tibble(
  x = rep(x, 10),
  true_y = intercept + slope * x,
  y = intercept + slope * x + rnorm(length(x), sd = 50),
  simulation = rep(1:10, each = length(x) / 10)
) 

animation_data %>% 
  ggplot(aes(x = x, y = y)) +
  geom_point() +
  geom_abline(slope = slope, intercept = intercept, col = "blue") +
  geom_smooth(method = "lm", se = FALSE, col = "red") +
  transition_states(states = simulation) +
  theme_bw()

slope_simulation <- function(
  regressor,
  intercept,
  slope,
  error_distribution = "normal",
  n_reps,
  alpha1,
  alpha2
) {
  n <- length(regressor) # Number of observations
  z <- (regressor - mean(regressor)) / sd(regressor) # Used for correlated errors
  count_in <- 0
  difference <- numeric(n_reps)
  for (i in 1:n_reps) {
    if (error_distribution == "normal") {
      error <- rnorm(n)
    } else if (error_distribution == "chisq") {
      error <- rchisq(n, 0.5) - 0.5
    } else if (error_distribution == "correlated") {
      error <- rnorm(n, 0.2 * z, 1)
    } else {
      stop("Unknown error distribution")
    }
    
    response <- intercept + slope * regressor + error
    fm <- lm(response ~ regressor)
    
    # Extract slope estimate
    fit_slope <- fm$coefficients[2]
    
    # Calculate difference between slope estimate and slope
    difference[i] <- fit_slope - slope
    
    # Confidence interval on fit_slope
    slope_conf <- fit_slope + c(-1, 1) * qt(1 - alpha1/2, fm$df.residual) * summary(fm)$coefficients[2, "Std. Error"]
    count_in <- count_in + (slope_conf[1] < slope && slope < slope_conf[2])
  }
  
  # Calculate coverage
  coverage <- count_in / n_reps
  
  # Confidence interval on coverage
  coverage_ci <- coverage + c(-1, 1) * qnorm(1 - alpha2/2) * sqrt((coverage  * (1 - coverage))/n_reps)
  
  # Bias estimate
  bias <- mean(difference)
  
  # Confidence interval on bias
  bias_ci <- bias + c(-1, 1) * qnorm(1 - alpha2/2) * sd(difference) / sqrt(n_reps)
  
  # Output
  list(
    coverage = c(
      lower = coverage_ci[1],
      estimate = coverage,
      upper = coverage_ci[2]
    ),
    bias = c(
      lower = bias_ci[1],
      estimate = bias,
      upper = bias_ci[2]
    )
  )
}

For each iteration of the n_reps iterations of the simulation, randomly generate the response (i.e., dependent) vector using the supplied intercept, slope, regressor, and error term distribution. Fit the simple linear regression model and compute a \(100(1-\alpha_1)\)% confidence interval on the slope parameter. Record whether it contains the true slope. Also, record the difference between the slope estimate and its true value.

The proportion of times that the confidence interval contains the true slope is a point estimate of the its coverage. Theory says that the coverage should be \(1-\alpha_1\) when the usual regression assumptions are met. In addition to providing a point estimate of the coverage, you will provide a \(100(1-\alpha_2)\)% confidence interval on the coverage. (Use the normal approximation to the binomial, which is justified by the Central Limit Theorem since nreps is large.)

The average difference between the slope estimate and its true value is a point estimate of the bias. In addition to providing a point estimate of the bias, you will provide a \(100(1-\alpha_2)\)% confidence interval on the bias. (Again, the Central Limit Theorem is applicable.)

The function should return the following six elements:

  1. A point estimate of the bias in estimating the slope
  2. The lower bound of a \(100(1-\alpha_2)\)% confidence interval for the bias
  3. The upper bound of a \(100(1-\alpha_2)\)% confidence interval for the bias
  4. A point estimate of the coverage of the \(100(1-\alpha_1)\)% confidence interval estimator of the slope
  5. The lower bound of a \(100(1-\alpha_2)\)% confidence interval for the coverage
  6. The upper bound of a \(100(1-\alpha_2)\)% confidence interval for the coverage

Evaluate

Consider the following questions:

simulation <- function(
  x, 
  intercept = 7, 
  slope = 0.5, 
  n_reps = 5000, 
  alpha1 = 0.10, 
  alpha2 = 0.05) {
  for (error_dist in c("normal", "chisq", "correlated")) {
    print(paste("Results for", error_dist, "errors:"))
    results <- slope_simulation(x, intercept, slope, error_dist, n_reps, alpha1, alpha2)
    print(paste0("Bias:     ", results$bias[2], " (", results$bias[1], ", ", results$bias[3], ")"))
    print(paste0("Coverage: ", results$coverage[2], " (", results$coverage[1], ", ", results$coverage[3], ")"))
  }
}

simulation(x = seq(-2, 2, length = 50))
[1] "Results for normal errors:"
[1] "Bias:     -0.000847634253262622 (-0.00418252079468957, 0.00248725228816433)"
[1] "Coverage: 0.9008 (0.89251422542325, 0.90908577457675)"
[1] "Results for chisq errors:"
[1] "Bias:     0.00181781368695764 (-0.00154552069960945, 0.00518114807352473)"
[1] "Coverage: 0.8946 (0.886088661926246, 0.903111338073754)"
[1] "Results for correlated errors:"
[1] "Bias:     0.167612312729427 (0.164251314366211, 0.170973311092642)"
[1] "Coverage: 0.5998 (0.586219840783989, 0.613380159216011)"
simulation(x = seq(-2, 2, length = 3))
[1] "Results for normal errors:"
[1] "Bias:     0.000966198337616723 (-0.00884107296014191, 0.0107734696353754)"
[1] "Coverage: 0.8968 (0.888367608976036, 0.905232391023964)"
[1] "Results for chisq errors:"
[1] "Bias:     -0.00725001709345427 (-0.0167822938467396, 0.00228225965983101)"
[1] "Coverage: 0.9486 (0.942479509777187, 0.954720490222813)"
[1] "Results for correlated errors:"
[1] "Bias:     0.103372396296765 (0.0934498034634754, 0.113294989130054)"
[1] "Coverage: 0.8942 (0.885674433319529, 0.902725566680471)"
LS0tCnRpdGxlOiAiU2ltdWxhdGlvbiBTdHVkeSBmb3IgU2ltcGxlIExpbmVhciBSZWdyZXNzaW9uIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpgYGB7ciBzZXR1cCwgaW5jbHVkZT1GQUxTRX0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KGVjaG8gPSBUUlVFLCBjYWNoZSA9IFRSVUUsIG1lc3NhZ2UgPSBGQUxTRSkKYGBgCgpSZWNhbGwgdGhhdCB0aGUgc2ltcGxlIGxpbmVhciByZWdyZXNzaW9uIG1vZGVsIGlzOgoKJCQKeV9pID0gXGJldGFfMCArIFxiZXRhXzEgeF9pICsgXGVwc2lsb25faSBccXVhZCBcdGV4dHtmb3J9IH4gaT0xLFxsZG90cyxuLAokJAoKd2hlcmUgdGhlIGVycm9yIHRlcm1zICRcZXBzaWxvbl8xLFxsZG90cyxcZXBzaWxvbl9uJCBhcmUgbWVhbiB6ZXJvLCBpbmRlcGVuZGVudCwgCmFuZCBpZGVudGljYWxseSBkaXN0cmlidXRlZCByYW5kb20gdmFyaWFibGVzLgoKTGV0ICRcaGF0e1xiZXRhfV8xJCBkZW5vdGUgdGhlIHBhcmFtZXRlciBlc3RpbWF0ZSBvZiAkXGJldGFfMSQuIFRoZSBiaWFzCiRcaGF0e1xiZXRhfV8xJCBpbiBlc3RpbWF0aW5nICRcYmV0YV8xJCBpcyBkZWZpbmVkIGFzIHRoZSBleHBlY3RlZCAoaS5lLiwgbWVhbikKdmFsdWUgb2YgJFxoYXR7XGJldGF9XzEgLSBcYmV0YV8xJC4gQSAkMTAwKDEtXGFscGhhXzEpJCUgY29uZmlkZW5jZSBpbnRlcnZhbCBmb3IKdGhlIHNsb3BlICRcYmV0YV8xJCAoYXNzdW1pbmcgaW5kZXBlbmRlbnQsIG5vcm1hbGx5IGRpc3RyaWJ1dGVkIGVycm9ycykgaXMKJFxoYXR7XGJldGF9XzEgXHBtIHRfe1xhbHBoYV8xLzIsXHRleHR7ZGZ9fSB+IFxoYXR7XHNpZ21hfV97XGhhdHtcYmV0YX1fMX0kLAp3aGVyZSAkXHRleHR7ZGZ9JCBpcyB0aGUgcmVzaWR1YWwgKGEuay5hLiwgZXJyb3IpIGRlZ3JlZXMgb2YgZnJlZWRvbSwKJFxoYXR7XHNpZ21hfV97XGhhdHtcYmV0YX1fMX0kIGlzIHRoZSBzdGFuZGFyZCBlcnJvciBvZiB0aGUgc2xvcGUgZXN0aW1hdGUsIGFuZAokdF97XGFscGhhXzEvMixkZn0kIGlzIHRoZSAkXGFscGhhXzEvMiQgdXBwZXIgcXVhbnRpbGUgb2YgdGhlICR0JC1kaXN0cmlidXRpb24Kd2l0aCAkXHRleHR7ZGZ9JCBkZWdyZWVzIG9mIGZyZWVkb20uIFRoZSBjb3ZlcmFnZSBvZiBhIGNvbmZpZGVuY2UgaW50ZXJ2YWwKZXN0aW1hdG9yIGlzIHRoZSBsb25nLXJ1biByZWxhdGl2ZSBmcmVxdWVuY3kgdGhhdCB0aGUgcHJvY2VkdXJlIGdlbmVyYXRlcwppbnRlcnZhbHMgdGhhdCBjb250YWluIHRoZSB0cnVlIHZhbHVlLgoKVGhlIHJlZ3Jlc3Npb24gY29lZmZpY2llbnRzICRcYmV0YV8wJCBhbmQgJFxiZXRhXzEkIGNhbiBiZSBlc3RpbWF0ZWQgdXNpbmcgdGhlCm1ldGhvZCBvZiBsZWFzdC1zcXVhcmVzLiBJbiBSLCBpZiB4IGFuZCB5IGFyZSBudW1lcmljIHZlY3RvcnMgb2YgdGhlIHNhbWUKbGVuZ3RoLCB0aGUgbGVhc3Qtc3F1YXJlcyBwYXJhbWV0ZXIgZXN0aW1hdGVzLCBzdGFuZGFyZCBlcnJvcnMsIGRlZ3JlZXMgb2YKZnJlZWRvbSwgZXRjLiBhcmUgb2J0YWluZWQgdXNpbmc6CgpgYGB7ciwgZXZhbD1GQUxTRX0KZm0gPC0gbG0oeSB+IHgpCmZtCnN1bW1hcnkoZm0pCmBgYAoKV2hpbGUgdGhlb3J5IHByb3ZpZGVzIHJlc3VsdHMgb24gdGhlIGJpYXMgYW5kIGNvdmVyYWdlIGluIHNpbXBsZSBsaW5lYXIKcmVncmVzc2lvbiBtb2RlbHMsIHRoZSBnb2FsIGlzIHRvIHBlcmZvcm0gYSBzaW11bGF0aW9uIHN0dWR5IHRvIGludmVzdGlnYXRlOgoKKyBUaGUgYmlhcyBpbiBlc3RpbWF0aW5nIHRoZSBzbG9wZSwgYW5kCisgVGhlIGNvdmVyYWdlIG9mIHRoZSAkMTAwKDEtXGFscGhhXzEpJCUgY29uZmlkZW5jZSBpbnRlcnZhbCBlc3RpbWF0b3IgZm9yIHRoZQpzbG9wZSBiYXNlZCBvbiB0aGUgdXN1YWwgbm9ybWFsaXR5IGFzc3VtcHRpb24uCgpUaGUgZm9sbG93aW5nIGZ1bmN0aW9uIGNvbmR1Y3RzIHN1Y2ggYSBzaW11bGF0aW9uIHN0dWR5LiBUaGUgYXJndW1lbnRzIHRvIHRoZQpmdW5jdGlvbiBhcmU6CgorIGByZWdyZXNzb3JgOiBJbmRlcGVuZGVudCB2YXJpYWJsZSBwYXNzZWQgYXMgYSBudW1lcmljIHZlY3RvciBvZiBhcmJpdHJhcnkKbGVuZ3RoLgorIGBpbnRlcmNlcHRgOiBUcnVlIGludGVyY2VwdCBwYXNzZWQgYXMgYSBudW1lcmljIHZlY3RvciBvZiBsZW5ndGggb25lLgorIGBzbG9wZWA6IFRydWUgc2xvcGUgcGFzc2VkIGFzIGEgbnVtZXJpYyB2ZWN0b3Igb2YgbGVuZ3RoIG9uZS4KKyBgZXJyb3JfZGlzdHJpYnV0aW9uYDogRGlzdHJpYnV0aW9uIG9mIHRoZSBlcnJvciB0ZXJtIHBhc3NlZCBhcyBhIGNoYXJhY3Rlcgp2ZWN0b3Igb2YgbGVuZ3RoIG9uZSwgZXF1YWxpbmcgb25lIG9mOgogICAgKyAibm9ybWFsIjogRXJyb3IgdGVybXMgYXJlIGluZGVwZW5kZW50IGFuZCBzdGFuZGFyZCBub3JtYWxseSBkaXN0cmlidXRlZC4KICAgICsgImNoaXNxIjogRXJyb3IgdGVybSBhcmUgaW5kZXBlbmRlbnRseSBjaGktc3F1YXJlZCBkaXN0cmlidXRlZCB3aXRoIDAuNQogICAgZGVncmVlcyBvZiBmcmVlZG9tIGFuZCBzaGlmdGVkIHRvIHRoZSBsZWZ0IGJ5IDAuNSAoc28gdGhhdCB0aGUgbWVhbiBvZiB0aGUKICAgIGVycm9yIHRlcm0gaXMgemVybykuICBUaGlzIGhhcyB0aGUgZWZmZWN0IG9mIHZpb2xhdGluZyB0aGUgbm9ybWFsaXR5CiAgICBhc3N1bXB0aW9uIG9mIHRoZSBlcnJvciB0ZXJtcyAkXGVwc2lsb24kJ3MuCiAgICArICJjb3JyZWxhdGVkIjogRXJyb3IgdGVybSBmb3Igb2JzZXJ2YXRpb24gaSBpcyBub3JtYWxseSBkaXN0cmlidXRlZCB3aXRoCiAgICBtZWFuICQwLjIgel9pJCBhbmQgdmFyaWFuY2UgMSwgd2hlcmUgJHpfaSA9ICh4X2ktXGJhcnt4fSkvcyQgaXMgdGhlCiAgICBzdGFuZGFyZGl6ZWQgdmFsdWUgb2YgcmVncmVzc29yIGZvciBvYnNlcnZhdGlvbiBpLiAgVGhpcyBoYXMgdGhlIGVmZmVjdCBvZgogICAgbWFraW5nIHRoZSBlcnJvciB0ZXJtcyAkXGVwc2lsb24kJ3MgY29ycmVsYXRlZCwgdGhlcmVieSB2aW9sYXRpbmcgdGhlCiAgICBpbmRlcGVuZGVuY2UgYXNzdW1wdGlvbi4KKyBgbl9yZXBzYDogTnVtYmVyIG9mIHJlcGxpY2F0aW9ucy4KKyBgYWxwaGExYDogRXF1YWxzIDEuMCBtaW51cyB0aGUgY29uZmlkZW5jZSBjb2VmZmljaWVudCBmb3IgdGhlIGNvbmZpZGVuY2UKaW50ZXJ2YWwgZXN0aW1hdG9yIG9mIHRoZSBzbG9wZS4gRm9yIGV4YW1wbGUsIGZvciA5NSUgY29uZmlkZW5jZSBpbnRlcnZhbHMgb24KdGhlIHNsb3BlLCBhbHBoYTEgPSAwLjA1LgorIGBhbHBoYTJgOiBFcXVhbHMgMS4wIG1pbnVzIHRoZSBjb25maWRlbmNlIGNvZWZmaWNpZW50IGZvciB0aGUgY29uZmlkZW5jZQppbnRlcnZhbHMgb2YgdGhlIGJpYXMgYW5kIGNvdmVyYWdlIHByb2JhYmlsaXRpZXMuIEZvciBleGFtcGxlLCBmb3IgYSA5MCUKY29uZmlkZW5jZSBpbnRlcnZhbCBvbiB0aGUgYmlhcywgYWxwaGEyID0gMC4xMC4KCmBgYHtyfQooeCA8LSBzZXEoLTIsIDIsIGxlbmd0aCA9IDUwKSkKaW50ZXJjZXB0IDwtIDU3CnNsb3BlIDwtIDExMgp0cnVlX3kgPC0gaW50ZXJjZXB0ICsgc2xvcGUgKiB4CgpwbG90KHgsIHRydWVfeSwgJ2wnKQpgYGAKCgpgYGB7cn0KZXJyb3IgPC0gcm5vcm0obGVuZ3RoKHgpLCBzZCA9IDEwKQpzaW11bGF0ZWRfeSA8LSB0cnVlX3kgKyBlcnJvcgpgYGAKCgpgYGB7cn0KcGxvdCh4LCB0cnVlX3ksIHR5cGUgPSAnbCcpCnBvaW50cyh4LCBzaW11bGF0ZWRfeSkKYGBgCgpgYGB7cn0KKGxtX291dCA8LSBsbShzaW11bGF0ZWRfeSB+IHgpKQpgYGAKCmBgYHtyfQooYmlhcyA8LSBsbV9vdXQkY29lZmZpY2llbnRzWzJdIC0gc2xvcGUpCmBgYAoKYGBge3J9CmFscGhhIDwtIC4xCmNpIDwtIGxtX291dCRjb2VmZmljaWVudHNbMl0gKyBjKC0xLCAxKSAqIHF0KDEgLSBhbHBoYSwgZGYgPSBsbV9vdXQkZGYucmVzaWR1YWwpICogc3VtbWFyeShsbV9vdXQpJGNvZWZmaWNpZW50c1syLCAiU3RkLiBFcnJvciJdCmMoCiAgbG93ZXIgPSBjaVsxXSwKICBlc3RpbWF0ZSA9IGxtX291dCRjb2VmZmljaWVudHNbMl0sCiAgdXBwZXIgPSBjaVsyXQopCmBgYAoKYGBge3IsIGV2YWwgPSBGQUxTRX0KbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkoZ2dhbmltYXRlKQoKYW5pbWF0aW9uX2RhdGEgPC0gdGliYmxlKAogIHggPSByZXAoeCwgMTApLAogIHRydWVfeSA9IGludGVyY2VwdCArIHNsb3BlICogeCwKICB5ID0gaW50ZXJjZXB0ICsgc2xvcGUgKiB4ICsgcm5vcm0obGVuZ3RoKHgpLCBzZCA9IDUwKSwKICBzaW11bGF0aW9uID0gcmVwKDE6MTAsIGVhY2ggPSBsZW5ndGgoeCkgLyAxMCkKKSAKCmFuaW1hdGlvbl9kYXRhICU+JSAKICBnZ3Bsb3QoYWVzKHggPSB4LCB5ID0geSkpICsKICBnZW9tX3BvaW50KCkgKwogIGdlb21fYWJsaW5lKHNsb3BlID0gc2xvcGUsIGludGVyY2VwdCA9IGludGVyY2VwdCwgY29sID0gImJsdWUiKSArCiAgZ2VvbV9zbW9vdGgobWV0aG9kID0gImxtIiwgc2UgPSBGQUxTRSwgY29sID0gInJlZCIpICsKICB0cmFuc2l0aW9uX3N0YXRlcyhzdGF0ZXMgPSBzaW11bGF0aW9uKSArCiAgdGhlbWVfYncoKQpgYGAKCiFbXShyZWdyZXNzaW9uX2FuaW1hdGlvbi5naWYpCgpgYGB7cn0Kc2xvcGVfc2ltdWxhdGlvbiA8LSBmdW5jdGlvbigKICByZWdyZXNzb3IsCiAgaW50ZXJjZXB0LAogIHNsb3BlLAogIGVycm9yX2Rpc3RyaWJ1dGlvbiA9ICJub3JtYWwiLAogIG5fcmVwcywKICBhbHBoYTEsCiAgYWxwaGEyCikgewogIG4gPC0gbGVuZ3RoKHJlZ3Jlc3NvcikgIyBOdW1iZXIgb2Ygb2JzZXJ2YXRpb25zCiAgeiA8LSAocmVncmVzc29yIC0gbWVhbihyZWdyZXNzb3IpKSAvIHNkKHJlZ3Jlc3NvcikgIyBVc2VkIGZvciBjb3JyZWxhdGVkIGVycm9ycwogIGNvdW50X2luIDwtIDAKICBkaWZmZXJlbmNlIDwtIG51bWVyaWMobl9yZXBzKQogIGZvciAoaSBpbiAxOm5fcmVwcykgewogICAgaWYgKGVycm9yX2Rpc3RyaWJ1dGlvbiA9PSAibm9ybWFsIikgewogICAgICBlcnJvciA8LSBybm9ybShuKQogICAgfSBlbHNlIGlmIChlcnJvcl9kaXN0cmlidXRpb24gPT0gImNoaXNxIikgewogICAgICBlcnJvciA8LSByY2hpc3EobiwgMC41KSAtIDAuNQogICAgfSBlbHNlIGlmIChlcnJvcl9kaXN0cmlidXRpb24gPT0gImNvcnJlbGF0ZWQiKSB7CiAgICAgIGVycm9yIDwtIHJub3JtKG4sIDAuMiAqIHosIDEpCiAgICB9IGVsc2UgewogICAgICBzdG9wKCJVbmtub3duIGVycm9yIGRpc3RyaWJ1dGlvbiIpCiAgICB9CiAgICAKICAgIHJlc3BvbnNlIDwtIGludGVyY2VwdCArIHNsb3BlICogcmVncmVzc29yICsgZXJyb3IKICAgIGZtIDwtIGxtKHJlc3BvbnNlIH4gcmVncmVzc29yKQogICAgCiAgICAjIEV4dHJhY3Qgc2xvcGUgZXN0aW1hdGUKICAgIGZpdF9zbG9wZSA8LSBmbSRjb2VmZmljaWVudHNbMl0KICAgIAogICAgIyBDYWxjdWxhdGUgZGlmZmVyZW5jZSBiZXR3ZWVuIHNsb3BlIGVzdGltYXRlIGFuZCBzbG9wZQogICAgZGlmZmVyZW5jZVtpXSA8LSBmaXRfc2xvcGUgLSBzbG9wZQogICAgCiAgICAjIENvbmZpZGVuY2UgaW50ZXJ2YWwgb24gZml0X3Nsb3BlCiAgICBzbG9wZV9jb25mIDwtIGZpdF9zbG9wZSArIGMoLTEsIDEpICogcXQoMSAtIGFscGhhMS8yLCBmbSRkZi5yZXNpZHVhbCkgKiBzdW1tYXJ5KGZtKSRjb2VmZmljaWVudHNbMiwgIlN0ZC4gRXJyb3IiXQogICAgY291bnRfaW4gPC0gY291bnRfaW4gKyAoc2xvcGVfY29uZlsxXSA8IHNsb3BlICYmIHNsb3BlIDwgc2xvcGVfY29uZlsyXSkKICB9CiAgCiAgIyBDYWxjdWxhdGUgY292ZXJhZ2UKICBjb3ZlcmFnZSA8LSBjb3VudF9pbiAvIG5fcmVwcwogIAogICMgQ29uZmlkZW5jZSBpbnRlcnZhbCBvbiBjb3ZlcmFnZQogIGNvdmVyYWdlX2NpIDwtIGNvdmVyYWdlICsgYygtMSwgMSkgKiBxbm9ybSgxIC0gYWxwaGEyLzIpICogc3FydCgoY292ZXJhZ2UgICogKDEgLSBjb3ZlcmFnZSkpL25fcmVwcykKICAKICAjIEJpYXMgZXN0aW1hdGUKICBiaWFzIDwtIG1lYW4oZGlmZmVyZW5jZSkKICAKICAjIENvbmZpZGVuY2UgaW50ZXJ2YWwgb24gYmlhcwogIGJpYXNfY2kgPC0gYmlhcyArIGMoLTEsIDEpICogcW5vcm0oMSAtIGFscGhhMi8yKSAqIHNkKGRpZmZlcmVuY2UpIC8gc3FydChuX3JlcHMpCiAgCiAgIyBPdXRwdXQKICBsaXN0KAogICAgY292ZXJhZ2UgPSBjKAogICAgICBsb3dlciA9IGNvdmVyYWdlX2NpWzFdLAogICAgICBlc3RpbWF0ZSA9IGNvdmVyYWdlLAogICAgICB1cHBlciA9IGNvdmVyYWdlX2NpWzJdCiAgICApLAogICAgYmlhcyA9IGMoCiAgICAgIGxvd2VyID0gYmlhc19jaVsxXSwKICAgICAgZXN0aW1hdGUgPSBiaWFzLAogICAgICB1cHBlciA9IGJpYXNfY2lbMl0KICAgICkKICApCn0KYGBgCgpGb3IgZWFjaCBpdGVyYXRpb24gb2YgdGhlIGBuX3JlcHNgIGl0ZXJhdGlvbnMgb2YgdGhlIHNpbXVsYXRpb24sIHJhbmRvbWx5CmdlbmVyYXRlIHRoZSByZXNwb25zZSAoaS5lLiwgZGVwZW5kZW50KSB2ZWN0b3IgdXNpbmcgdGhlIHN1cHBsaWVkIGludGVyY2VwdCwKc2xvcGUsIHJlZ3Jlc3NvciwgYW5kIGVycm9yIHRlcm0gZGlzdHJpYnV0aW9uLiBGaXQgdGhlIHNpbXBsZSBsaW5lYXIgcmVncmVzc2lvbgptb2RlbCBhbmQgY29tcHV0ZSBhICQxMDAoMS1cYWxwaGFfMSkkJSBjb25maWRlbmNlIGludGVydmFsIG9uIHRoZSBzbG9wZQpwYXJhbWV0ZXIuIFJlY29yZCB3aGV0aGVyIGl0IGNvbnRhaW5zIHRoZSB0cnVlIHNsb3BlLiBBbHNvLCByZWNvcmQgdGhlCmRpZmZlcmVuY2UgYmV0d2VlbiB0aGUgc2xvcGUgZXN0aW1hdGUgYW5kIGl0cyB0cnVlIHZhbHVlLgoKVGhlIHByb3BvcnRpb24gb2YgdGltZXMgdGhhdCB0aGUgY29uZmlkZW5jZSBpbnRlcnZhbCBjb250YWlucyB0aGUgdHJ1ZSBzbG9wZSBpcwphIHBvaW50IGVzdGltYXRlIG9mIHRoZSBpdHMgY292ZXJhZ2UuIFRoZW9yeSBzYXlzIHRoYXQgdGhlIGNvdmVyYWdlIHNob3VsZCBiZQokMS1cYWxwaGFfMSQgd2hlbiB0aGUgdXN1YWwgcmVncmVzc2lvbiBhc3N1bXB0aW9ucyBhcmUgbWV0LiBJbiBhZGRpdGlvbiB0bwpwcm92aWRpbmcgYSBwb2ludCBlc3RpbWF0ZSBvZiB0aGUgY292ZXJhZ2UsIHlvdSB3aWxsIHByb3ZpZGUgYQokMTAwKDEtXGFscGhhXzIpJCUgY29uZmlkZW5jZSBpbnRlcnZhbCBvbiB0aGUgY292ZXJhZ2UuIChVc2UgdGhlIG5vcm1hbAphcHByb3hpbWF0aW9uIHRvIHRoZSBiaW5vbWlhbCwgd2hpY2ggaXMganVzdGlmaWVkIGJ5IHRoZSBDZW50cmFsIExpbWl0IFRoZW9yZW0Kc2luY2UgbnJlcHMgaXMgbGFyZ2UuKQoKVGhlIGF2ZXJhZ2UgZGlmZmVyZW5jZSBiZXR3ZWVuIHRoZSBzbG9wZSBlc3RpbWF0ZSBhbmQgaXRzIHRydWUgdmFsdWUgaXMgYSBwb2ludAplc3RpbWF0ZSBvZiB0aGUgYmlhcy4gSW4gYWRkaXRpb24gdG8gcHJvdmlkaW5nIGEgcG9pbnQgZXN0aW1hdGUgb2YgdGhlIGJpYXMsIHlvdQp3aWxsIHByb3ZpZGUgYSAkMTAwKDEtXGFscGhhXzIpJCUgY29uZmlkZW5jZSBpbnRlcnZhbCBvbiB0aGUgYmlhcy4gKEFnYWluLCB0aGUKQ2VudHJhbCBMaW1pdCBUaGVvcmVtIGlzIGFwcGxpY2FibGUuKQoKVGhlIGZ1bmN0aW9uIHNob3VsZCByZXR1cm4gdGhlIGZvbGxvd2luZyBzaXggZWxlbWVudHM6CgoxLiBBIHBvaW50IGVzdGltYXRlIG9mIHRoZSBiaWFzIGluIGVzdGltYXRpbmcgdGhlIHNsb3BlCjIuIFRoZSBsb3dlciBib3VuZCBvZiBhICQxMDAoMS1cYWxwaGFfMikkJSBjb25maWRlbmNlIGludGVydmFsIGZvciB0aGUgYmlhcwozLiBUaGUgdXBwZXIgYm91bmQgb2YgYSAkMTAwKDEtXGFscGhhXzIpJCUgY29uZmlkZW5jZSBpbnRlcnZhbCBmb3IgdGhlIGJpYXMKNC4gQSBwb2ludCBlc3RpbWF0ZSBvZiB0aGUgY292ZXJhZ2Ugb2YgdGhlICQxMDAoMS1cYWxwaGFfMSkkJSBjb25maWRlbmNlCmludGVydmFsIGVzdGltYXRvciBvZiB0aGUgc2xvcGUKNS4gVGhlIGxvd2VyIGJvdW5kIG9mIGEgJDEwMCgxLVxhbHBoYV8yKSQlIGNvbmZpZGVuY2UgaW50ZXJ2YWwgZm9yIHRoZSBjb3ZlcmFnZQo2LiBUaGUgdXBwZXIgYm91bmQgb2YgYSAkMTAwKDEtXGFscGhhXzIpJCUgY29uZmlkZW5jZSBpbnRlcnZhbCBmb3IgdGhlIGNvdmVyYWdlCgojIyBFdmFsdWF0ZQpDb25zaWRlciB0aGUgZm9sbG93aW5nIHF1ZXN0aW9uczoKCisgVW5kZXIgd2hhdCBlcnJvciB0ZXJtcyBpcyB0aGUgZXN0aW1hdG9yIG9mIHRoZSBzbG9wZSBiaWFzZWQ/CisgRG9lcyB0aGUgY29uZmlkZW5jZSBpbnRlcnZhbCBlc3RpbWF0b3IgaGF2ZSB0aGUgcmlnaHQgY292ZXJhZ2UgdW5kZXIgdGhlCm5vcm1hbGx5IGRpc3RyaWJ1dGVkIGVycm9yIHRlcm1zPworIFdoYXQgaXMgdGhlIGNvdmVyYWdlIHVuZGVyIHRoZSBjaGktc3F1YXJlZCBkaXN0cmlidXRlZCBlcnJvciB0ZXJtcz8KKyBXaGF0IGlzIHRoZSBjb3ZlcmFnZSB1bmRlciB0aGUgY29ycmVsYXRlZCBlcnJvciB0ZXJtcz8KKyBGb3IgdGhlIGVycm9yIHRlcm1zIGluIHdoaWNoIHRoZSBjb3ZlcmFnZSBpcyBub3QgY29ycmVjdCwgdW5kZXIgd2hhdCBzaXR1YXRpb24KaXMgaXQgbm90aWNlYWJsZT8gV2hlbiwgaWYgZXZlciwgZG9lcyB0aGlzIGNvdmVyYWdlIHByb2JsZW0gZ28gYXdheT8gV2hhdApwaGVub21lbm9uIG1ha2VzIGl0IGdvIGF3YXk/CgpgYGB7cn0Kc2ltdWxhdGlvbiA8LSBmdW5jdGlvbigKICB4LCAKICBpbnRlcmNlcHQgPSA3LCAKICBzbG9wZSA9IDAuNSwgCiAgbl9yZXBzID0gNTAwMCwgCiAgYWxwaGExID0gMC4xMCwgCiAgYWxwaGEyID0gMC4wNSkgewogIGZvciAoZXJyb3JfZGlzdCBpbiBjKCJub3JtYWwiLCAiY2hpc3EiLCAiY29ycmVsYXRlZCIpKSB7CiAgICBwcmludChwYXN0ZSgiUmVzdWx0cyBmb3IiLCBlcnJvcl9kaXN0LCAiZXJyb3JzOiIpKQogICAgcmVzdWx0cyA8LSBzbG9wZV9zaW11bGF0aW9uKHgsIGludGVyY2VwdCwgc2xvcGUsIGVycm9yX2Rpc3QsIG5fcmVwcywgYWxwaGExLCBhbHBoYTIpCiAgICBwcmludChwYXN0ZTAoIkJpYXM6ICAgICAiLCByZXN1bHRzJGJpYXNbMl0sICIgKCIsIHJlc3VsdHMkYmlhc1sxXSwgIiwgIiwgcmVzdWx0cyRiaWFzWzNdLCAiKSIpKQogICAgcHJpbnQocGFzdGUwKCJDb3ZlcmFnZTogIiwgcmVzdWx0cyRjb3ZlcmFnZVsyXSwgIiAoIiwgcmVzdWx0cyRjb3ZlcmFnZVsxXSwgIiwgIiwgcmVzdWx0cyRjb3ZlcmFnZVszXSwgIikiKSkKICB9Cn0KCnNpbXVsYXRpb24oeCA9IHNlcSgtMiwgMiwgbGVuZ3RoID0gNTApKQpgYGAKCmBgYHtyfQpzaW11bGF0aW9uKHggPSBzZXEoLTIsIDIsIGxlbmd0aCA9IDMpKQpgYGAKCg==